Anisotropic spatially heterogeneous dynamics on the a and (3 relaxation time scales 

studied via a four-point correlation function 
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We examine the anisotropy of a four-point correlation function G4,{k,f;t) and it's associated 
structure factor S4,{k,q\t) calculated using Brownian Dynamics computer simulations of a model 
glass forming system. These correlation functions measure the spatial correlations of the relaxation 
of different particles, and we examine the time and temperature dependence of the anisotropy. We 
find that the anisotropy is strongest at nearest neighbor distances at time scales corresponding to 
the peak of the non-Gaussian parameter oi2{t) — 3(5r^(f))/[5((5r'^(t))'^] — 1, but is still pronounced 
around the a relaxation time. We find that the structure factor S4,{k,q;t) is anisotropic even for 
the smallest wave vector accessible in our simulation suggesting that our system (and other systems 
commonly used in computer simulations) may be too small to extract the g — > limit of the structure 
factor. We find that the determination of a dynamic correlation length from 5*4 (fc, q; t) is influenced 
by the anisotropy. We extract an effective anisotropic dynamic correlation length from the small q 
behavior of Si{k,q;t). 

PACS numbers: 



I. INTRODUCTION 



It is now generally aceepted that upon approaching the 
glass transition, the liquid's dynamics are becoming in- 
creasingly heterogeneous P, 0, However, the details 
of the spatial and temporal characteristics of dynamic 
heterogeneities are still being debated. In particular, 
the connection between heterogeneous dynamics and a 
growing dynamic correlation length has been the topic 
of many simulations B, j ^, B B fiol fill and a few 
experimental studies |12l . Il3l . Il4| . Four-point correlation 
functions have been introduced to facilitate the quantita- 
tive description of heterogeneous dynamics. The analy- 
sis of the spatial decay of these correlation functions was 
used to extract a dynamic correlation length. Recently, 
the mode-coupling theory has been extended and a the- 
oretical treatment of four- point correlation functions is 
starting to emerge 0, [H, [H, [H, [13, [III- However, in 
most simulation studies and in some theoretical treat- 
ments these four-point correlation functions have been 
assumed to be isotropic or they are isotropic by design. 

Researchers have noticed anisotropy in the correlated 
motion of particles on the (3 relaxation time scale, and re- 
cently this anisotropic motion has also been reported on 
the a relaxation time scale [l^. Doliwa and Heuer [2^ 
reported anisotropic correlated motion in a hard sphere 
system on the (3 relaxation time scale. Anisotropic mo- 
tion has also been extensively studied by Donati et al. 
and Gebremichael et al. [1, [2l[ who described the motion 
of "mobile" particles as "string-like", with mobile par- 
ticles following each other in one dimensional " strings" . 
Weeks et al. j22l | have reported anisotropic dynamics as- 
sociated with the break down of the "cage" surrounding 
a particle. They found that the correlations of the par- 
ticle's displacements depends on the initial separation of 
the particles. While particles that start at a separation 



corresponding the the first peak of the pair correlation 
function are most likely to move in the same direction, 
particles that start at a separation corresponding to the 
first minimum are more likely to initially move in oppo- 
site directions. 



In view of the experimental and simulational evidence 
for anisotropic correlations of particle's displacements, 
it should not be a surprise that four-point correlation 
functions designed to study these dynamics can also be 
anisotropic. However, this anisotropy is normally stud- 
ied for times less than the a relaxation time, thus it is 
uncertain if understanding this anisotropy is important 
for the structural relaxation of the liquid. Previously [l^ 
we reported on a four-point correlation function that is 
anisotropic on the a relaxation time scale as well as the 
(3 relaxation time scale for a model glass forming liq- 
uid. Since the spatial decay of this correlation function 
can be used to determine a dynamic length scale, the 
anisotropy introduces a complication in determining this 
length scale. 



In this paper we expand on previous work [19| . Af- 
ter describing the simulation in Sec. |TT1 we explore the 
anisotropic correlated dynamics by examining a four- 
point correlation function Gi{k,f]t), Sec. IIIIl and the 
associated structure factor Si{k^ q; t), Sec. lIVI We exam- 
ine the anisotropy at around nearest neighbor distances, 
which corresponds to local rearrangement of particles 
and its cage, and at large distances. We examine how 
the anisotropy influences the determination of a growing 
length scale accompanying the glass transition, and de- 
termine an effective anisotropic correlation length. We 
finish with a discussion of the results in Sec. IVl 
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II. SIMULATION 



relaxation of different particles. Consider the function 



We performed Brownian dynamics simulations of an 
80:20 binary mixture of 1000 particles introduced by 
Kob and Andersen [H, [l^. The interaction potential 
is Vapir) = 4eQ^[(cra/3/r)^^ - {aap/r)*^] where a, (3 G 
{A,B}, eAA = 1-0, Cab = 1-5, ess = 0.5, aAA = 1-0, 
o'AB = 0.8, and obb ~ 0.88 and the interaction poten- 
tial is cut at 2.5 Uap- Periodic boundary conditions were 
used with a box length of 9.4 uaa- The equation of mo- 
tion for the position of particle i is 



F,{t)+fj,{t), 



(1) 



where = 1-0 is the friction coefficient of an isolated 
particle and the force acting on a particle i is 



(2) 



with Vi being the gradient operator with respect to r^. 
The random force f](t) satisfies the fluctuation dissipation 
relation 



(3) 



where Do = ksT/^Q, is Boltzmann's constant, and 1 
is the unit tensor. The results are presented in terms of 
reduced units with aAA, ^aaI^b, and (j\a^oI^aa being 
the units of length, energy, and time, respectively. Since 
the equation of motion allows for diffusion of the center 
of mass, all results are presented relative to the center of 
mass. 

We present results for temperatures T = 0.45, 0.47, 
0.5, 0.55, 0.6, 0.8, 0.9, and 1.0. The onset of supercool- 
ing is around T = 1.0 and we use — 0.435 as the 
mode coupling temperature. As a means to expand the 
temperature scale, we will plot some quantities versus 
e = (T — T(.)/Tc. The equation of motion was integrated 
using a Heun algorithm with a small time step of 5 x 10~^. 
We ran an equilibration run that was at least half as 
long as a production run, and four production runs at 
each temperature. The results are an average over the 
production runs. We present results only for the larger 
and more abundant A particles. We define the a relax- 
ation time Tq, as through relation Fs{k]Ta) ~ e^^ for a 
wave vector around the first peak of the partial static 
structure factor for the A particles, which corresponds to 
|fc| = 7.25. 



III. FOUR-POINT CORRELATION FUNCTION 

G4{k,r-t) 

A. Definition and connection with overlap 
correlations 

We study a four-point correlation function that mea- 
sures the spatial and temporal correlations between the 



(4) 



where f„(i) is the position of particle n at a time t. The 
ensemble average of Fn{k; t) is the self- intermediate scat- 
tering function Fs{k]t), thus we will term Fn{k;t) the 
microscopic self-intermediate scattering function. The 
four-point correlation function 

V 



Gi{k,r-t) 



N2 



{Fn{k;t)Fmi-k;t)S[f-r'nmm) 



(5) 

measures the correlations between the microscopic self- 
intermediate scattering function at time t, pertaining to 
particles that are separated by a vector r at the initial 
time. In Eq. (O Turn — Tn ~ Tim V IS the volumc, and 
N is the number of particles. Notice that Gi{k,r]Q) = 
g{r) where _g(r) is the pair correlation function. In this 
work we choose |fc| to have the same value as the one 
that determines the a relaxation time, i.e. \k\ is located 
around the first peak of the partial static structure factor 
for the A particles. |fc | = 7.25. 

It should be noted that Gi{k,r;t) is, in general, com- 
plex. Its real and imaginary parts can be written in the 
following form 



Re[Gi{k,r\t)] = 
cosjfc 



(6) 



-y 



iV2 ^ 

Im[Gi{k,r;t)] 



[t) - rnm{0)]}S['r - r„ 



(7) 



7V2 



sin{fc ■ 



mm 



Eqs. ([7][8]) show that particles which are getting closer 
together or farther apart along the direction of vector 
k (i.e. are moving in the opposite direction or in the 
same direction along k) make the same contribution to 
the real part of Gi{k,r]t) but opposite contributions to 
its imaginary part. In particular, particles moving farther 
apart along the direction of vector k make a negative 
contribution to the imaginary part of G4(k,r]t). 

In several other simulational and experimental studies 
m, [IHl four-point correlation functions involving 
single-particle overlaps rather than the microscopic self- 
intermediate scattering functions were investigated. For 
example, Lacevic et. al Q used the following function 

M 

gf{r-t) = ^Y. Ma;t)wra{a;t)6[r~rnram). (8) 

where it;„(a; t) is the overlap function pertaining to par- 
ticle n. 



(9) 
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We would like to point out that g4{r; t) can be expressed 
in terms of functions which are generalizations of our 
G4{k,r;t), 

(10) 

where Gi^ki, k2, r; t) is defined as the correlation function 
of the microscopic self- intermediate scattering function at 
time t and calculated for different wave vectors, 

QA{kiM,r-t) = ^ ^ {Fn{ki]t)F,^{k2]t)5[r-rnrrM]), 

n^rn 

(11) 

and /(fc; a) = ATra'^ji{ka)/k with ji denoting a spherical 
Bessel function of the first kind. 

The present work is mostly concerned with the 
anisotropic nature of dynamic heterogeneities, which can 
be monitored using the four-point correlation function 
given by Eq. ((S]). In this context we would like to em- 
phasize that in principle the more general function (jlip 
is also anisotropic. However, any trace of this anisotropy 
is lost after the integration over wave vectors fci and ^2 
and thus the overlap correlation function ^ is, by con- 
struction, isotropic. 

B. Anisotropy of Gi{k,f\t) 

Since the functions F„{k;t) are sensitive to displace- 
ments of particles along the direction of k, then G4{k, r; t) 
measures interparticle correlations weighted by the dis- 
placements along the vector k. Particles which move in 
the direction perpendicular to k make a contribution to 
G4{k, r; t) which is the same as their contribution to the 
pair correlation function g{r). We notice that for t > 
four-point function Gi{k,r]t) is not isotropic, but de- 
pends on the angle 9 between r and k. Shown in the up- 
per figure in Fig.[T]is the real part 6*4 (fc, r; t) for T = 0.45 
calculated at t = Tq, and the lower figure shows the 
imaginary part. The maximum value of the real part 
of Gi{k,r;Ta) occurs for values of cos{9) corresponding 
to 6 ~ Q° and 9 — 180°, which shows that the correla- 
tions are most pronounced for f parallel and antiparallel 
to k. Thus, the correlations of the microscopic relaxation 
function is anisotropic on the the a relaxation time scale 
and the correlations are the strongest when neighboring 
particles move in the same or in opposite directions. 

To examine these anisotropic correlations at length 
scales around nearest neighbor distances, we expand 
G4(fc,r;t) into the Lcgendre polynomials 

G4(fc,f;t) =^L„(fc,r;0/'„(^^, (12) 

n 

where P„ is the nth Lcgendre Polynomial, k — k/k,f — 




FIG. 1: The real part of the correlation function G4{k,f;Ta) 
(upper figure) and the imaginary part of G4{k,f;Ta) (lower 
figure) for T=0.45 calculated at the a relaxation time. 

r/r, and 

Ln{k,r;t) = I Gi{k,r-t)P^{k-T)dr. (13) 

If Gi{k, f; t) docs not depend on the angle between k and 
r, then Ln{k, r; t) is zero for all n not equal to zero. Since 
there are nonzero real and imaginary parts to G^^k, r; t) 
for t > 0, then there are nonzero real and imaginary parts 
to L„(fc, r;i). By symmetry, the imaginary part is zero 
for even n, and the real part is zero for odd n. 

Shown in Fig. [2] is the real part Ln{k, r; r^) for n — 
and 2, and the imaginary part for n = 1 at the alpha 
relaxation time Tq, for T = 0.45. There is a peak in 
L2{k,r;Ta) and LQ{k,r;Ta) around the first peak of the 
pair correlation function g{r). The dashed lines in the 
figure are g(r)e~^. Note that due to our definition of 
the a relaxation time is the asymptotic limit of the 
isotropic component Lq a.t t = Ta, limj-^oo LQ{k,r;Ta) = 
Fg{k,Ta) = e~^. The positive peak in L2{k,r;Ta) in- 
dicates that particles that are initially separated by a 
distance corresponding to the first peak of g{r) have a 
tendency to move in the same direction or in opposite 
directions, while the values close to zero around the first 
minimum of the static structure factor can result from 
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FIG. 2: The real part of Ln{k, r; Tq,) for n = and 2, and the 
imaginary for n = 1 for T=0.45 calculated at the a relaxation 
time. The dashed line in the figures is g{r)e~^ where g{r) is 
the pair correlation function. 



motion which is perpendicular to the initial separation 
vector. The spatial variation of the correlated motion on 
these length scales has been reported previously in col- 
loidal suspensions [l^ and is related to the break up of 
the cage surrounding a particle. 

The variation of the imaginary part of Li(k,r]Ta) in- 
dicates that particles closer than the first peak of g{r) are 
more likely to move apart, while particles at a distance 
greater than this peak are more likely to move closer to- 
gether. In general, negative values of ii(fc, r; t) indicates 
that particles move farther apart while positive values 
indicate that particles move closer together. 

To look at the time dependence of the anisotropy, 
we calculated the height of the first peak of L2{k,r;t) 
as a function of time, which is shown in Fig. [3] for 
T = 1.0, 0.9, 0.8, 0.6, 0.55, 0.5, 0.47 and 0.45. The 
peak height starts at zero since the liquid is isotropic, 
then increases, reaches a maximum, and finally decreases 
to zero at long times. The height of the peak, tl2, 
is around the a relaxation time for high temperatures, 
Fig. m but its position increases slower with decreas- 
ing temperature than the a relaxation time and approx- 
imately follows the temperature dependence of the time 




FIG. 3: The time dependence of the the first peak of 
L2{k,r;Tc) for T = 1.0, 0.9, 0.8, 0.6, 0.55, 0.5, 0.47, 0.45, 
shown from left to right. 




(T - TJfl^ 



FIG. 4: The time at which the first peak of 1/2 (fe, r; t) reaches 
its maximum value, tl2, (squares) compared to the a relax- 
ation time, Ta, (circles) and the peak time of the standard 
non-gaussian parameter, Tng, (triangles). 



corresponding to the peak position of the standard non- 
Gaussian parameter a2{t) = 3{dr'^{t))/[5{5r'^{t))^] - 1, 
Tng (triangles in Fig. |4|). Furthermore, the maximum 
value does not monotonically increase with a decrease in 
the temperature, but rather reaches a maximum around 
T = 0.55, then begins to decrease with decreasing tem- 
perature. Thus the anisotropy around nearest neighbor 
distances initially increases upon supercooling the liquid, 
but reaches a maximum and begins to slowly decrease 
when the liquid is cooled further. It is not known if 
the peak height continues to decrease or saturates at low 
temperatures. 
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FIG. 5: The four-point correlation function S4,{k,q;t) for 6 = 
0, where is the angle between k and q, calculated at t = 0, 
0.1 Tq, Tc and IOtq for a temperature T = 0.45. 



FIG. 6: The four-point correlation function S4{k,q;Ta) for 
9 — 0° , 90°, and 180°, where 9 is the angle between q and k 
calculated for T = 0.45. 



IV. FOUR-POINT STRUCTURE FACTOR 

S4k,q;t) 

A. Anisotropy of S4(fc, q; t) 

To investigate the correlations between microscopic 
self-intermediate scattering functions at larger distances, 
we examined the the structure factor corresponding to 
G4(fc,r;i), 

N ^ 

Si{k,q-t) = l + -Hi{k,q-t) (14) 

where Hi{k, q; t) is the Fourier transform of G^ik, f; t) — 
F^{k-t). For q:^ 

S^ik,q;t) = l^(/'„(fc;t)/',„(-fc;t)e-«-™(0)). (15) 

Again, we fix |fc| to be around the position of the first 
peak of the static structure factor for the A particles, 
=7.25. 

Functions similar to p5p have been used to examine 
a growing dynamic length scale in glass forming liquids 
0, H, [13, El, In Fig. O we show results similar to 
those presented in, e.g. Ref. Specifically, we show in 
Silk, q^\;t) for T = 0.45 at times t ~ 0, O.Itq,, Tq, and 

lOTa. Note that for t = 0, Siik, q; 0) = S{q) where S{q) is 
the static structure factor for the A particles. We would 
like to emphasize that results shown in Fig. [5] are for one 
specific angle between k and q; the angle between q and 
k is zero. It should be noted that for this angle between 
vectors q and k, S4 does not depend on time for q = k. 
This follows from definition Si{k, k; t) — S{k) at all 
times. 



The usual interpretation of results shown in Fig. [5] is 
that the increase of S^i^k, q; t) at small q values suggests a 
growing dynamic length scale £,{t). To find the dynamic 
length scale, it is common to fit the small q behavior to a 
functional form and to examine the scaling of S4{k,q;t) 
for small q. In such a procedure it is implicitly assumed 
that S4{k,q]t) is isotropic. 

However, we find that S4{k,q;t) is not isotropic and 
depends on the angle between k and q. Shown in Fig. [5] 
is S4{k, g, Ta) for T = 0.45 and for 6 = 0° , 90°, and 180° 
where 9 is the angle between k and q. The anisotropy of 
Silk, q; t) adds a complication in finding a unique (,{t). 

Since we do not expect any slowly-decaying with 
increasing distance spatial correlations between self- 
intermediate scattering functions pertaining to different 
particles, we can safely assume that the (f — > limit 
of S/^{k,q;t) is well defined and it does not depend on 
the angle between vectors q and k. However, the re- 
sults shown in Fig. [5] suggest that the correlation length 
may be anisotropic. We would like to emphasize that 
our results are consistent with such a possibility but do 
not prove it. To prove that the correlation length is 
anisotropic one would need to simulate bigger systems 
in order to be able to examine the structure factor at 
smaller wave vectors q. 

We examine the anisotropy of the four-point structure 
factor by calculating the projection of S/^{k, q; t) onto the 
Legendre polynomials, 

2n + 1 P 

/„(fc, q; t) = ^— / Siik, q; t)Pn{k ■ q}dq. (16) 

Shown in Fig.[7{a) is /o(fc, q] Ta) (i.e., the angular average 
of S4) for T = 1.0, 0.8, 0.6, 0.55 and 0.45. In most sim- 
ulational studies of four-point correlation functions the 
correlation functions are shown as averages over different 
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FIG. 7: The projections Io{k,q;t), (a), and l2{k,q;t), (b), as 
described in the text. The projection /o(A;, q; t) is the average 
over angles 6 between k and q of S4,{k,q;t). If S4,{k,q;t) does 
not depend on 0, then l2{k,q;t) would be zero. Shown in 
the inset is l2{k,qo;ra) where go is the smallest wave vector 
allowed due to periodic boundary conditions as a function of 
temperature. The symbols in (a) and (b) correspond to the 
same temperatures. 



directions of wave vector q, thus the resuhs are similar 
to what is shown in Fig. \T[^a.). Note, however, that an 
average over different directions of q may not correspond 
to an angular average if the same number of wave vectors 
corresponding to each angle between q and k are not used 
in the average. Therefore, different routines to determine 
S4^{k,q;t) can lead to different conclusions, and our re- 
sults demonstrate that the averaging procedure needs to 
be performed with caution. 

Shown in Fig. [TJb) is l2{k,q;Ta) for T = 1.0, 0.8, 0.6, 
0.55 and 0.45. The non-zero values of I2 is a consequence 
of S4{k) being anisotropic on the a relaxation time scale. 
The anisotropy is largest for the smallest q values. The 
temperature dependence of I2 {k, qo]Ta) is shown as an in- 
set to Fig.[7]Jb). The anisotropy at the a relaxation time 
for qo grows with decreasing temperature until around 



FIG. 8: Time dependence of Si{k,qo;t) for different angles 
between k and go calculated for T — 0.45. The solid lines 
correspond to 61 = 0°, 30°, 45°, 60°, and 90° listed from 
bottom to top. The dashed lines are 120° , 135° , 150° and 180° 
listed from top to bottom. The vertical dotted line indicates 
the a relaxation time. 




FIG. 9: Time dependence of l2{k, go; t) where go is the small- 
est wave vector allowed due to periodic boundary conditions 
for T = 1.0, 0.9, 0.8, 0.6, 0.55, 0.5, 0.47, and 0.45 listed from 
left to right. 



T = 0.5, then it remains approximately constant. 



B. Time dependence of the anisotropy of S4{k,qo;t) 

We now turn to the examination of the time depen- 
dence of the anisotropy of Si{k, q; t). To this end, we set 
\q\ to be equal to the smallest wave vector allowed for our 
finite size simulation box, \q\ = qo = Stt/L, and calculate 
Si{k, qQ]t) as a function of time for different angles be- 
tween k and q. Results for T = 0.45 are shown in Fig. [51 
and the vertical line marks the a relaxation time. We see 
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FIG. 10: The time corresponding to the maximum value ol 
the magnitude of l2{k,qo;t), r/2, (circles) compared to the 
peak position of the non-Gaussian parameter T„g (diamonds) 
and the a relaxation time (squares). 



that Si{k,q]t) grows with increasing time, then reaches 
a maximum that depends on 9 for a time around the a 
relaxation time and finally decays to one at long times. 
Note that, while the position of the maximum is around 
the a relaxation time, the specific time at which the peak 
is reached depends on the angle between k and q. 

To determine the time dependence of the anisotropy, 
we examined l2{k,qo]t) where go is the smallest wave 
vector allowed due to periodic boundary condition, qa = 
2tt/L. As seen in Fig. [31 l2{k,qo;t) is zero at short and 
long times, but develops a peak at intermediate times. 
Note that the shape of l2{k,qo;t) is somewhat similar 
to that of L2{k,rpeak'it) shown in Fig. [3] except that 
l2{k,qo;t) is negative (the last fact could be expected 
from the relation between L2{k, r; t) and hik, q; t)). The 
peak height increases with decreasing temperature until 
T = 0.47, where it starts to decrease. However, as we 
show in the next subsection, the correlation length ob- 
tained from the fits at T = 0.45 are all close to or greater 
than half the box length, and it is currently unknown if 
the decrease in the peak height is a finite size effect. 

To determine when the anisotropy is a maximum 
at large distances, we found the time when l2{k,qo;t) 
reaches its maximum value, r/2. Shown in Fig. [TU] is 
the temperature dependence of r/2 (circles) compared to 
Ta (squares) and the peak position of the standard non- 
Gaussian parameter r„g (diamonds). We notice similar 
trends as with the time corresponding to the maximum 
value of L2{k, Vmax] t) in that the r/2 occurs around r„g 
and has a similar temperature dependence. 



C. Effective dynamic correlation length 

There has been some effort to determine the dy- 
namic correlation length by fitting functions similar to 



S4{k,q;t) to different functional forms [1, [iT|. Lacevic et 
al. used an Ornstein-Zernicke form A/ {I + {^q)^) to 
fit an overlap function Sf{q) that is isotropic by design, 
while Toninelh et al. [U used {A - C)/(l + (fg)'^) + C 
to fit a function similar to the one studied in this work. 
Lacevic et al. found a correlation length growing with 
time until the peak time in the associated four-point sus- 
ceptibility, and then decreasing. In contrast, Toninelli 
et al. found a correlation length growing with time even 
after the peak in the associated susceptibility. It is possi- 
ble that the difference between these findings was related 
th the presence of the new parameter (3 in the fit used 
by Toninelli et al. More recently, Berthier et al. @ used 
A/{1 + {£,q)^) and found a value oi (3 ~ 2.4 provided 
good fits to the same correlation function studied in the 
Ref. Q. Here we focus on a possible anisotropy of the 
correlation length at the time equal to the a relaxation 
time and we leave its time dependence for a future study. 
We started with 



Si{k,q]Ta) 



Si{k,Q-T^)-C 
1 + {^eq? + iaq)* 



(17) 



as a fitting function to extract the dynamic correlation 
length . In Eq. ([TT]) we added a constant C because of 
the growing baseline which can be seen in Fig.O We note 
that since we do not expect any slowly decaying spatial 
correlations, in the limit q Q, Si{k, q; Ta) should be in- 
dependent on the angle between k and q. In contrast, in 
Eq. p?|) we allowed for the dependence of the dynamic 
correlation length on the angle 6 between k and q. 
While fits to Eq. [T7] were very good for g < 3, the re- 
sults were not satisfactory. The values of S4{k, 0; t) were 
not consistent for different angles 9 between k and q and 
the length scales £,g obtained from the fits were greater 
than 40 at the lowest temperatures. To solve these prob- 
lems we performed the procedure described below. We 
emphasize that simulations of larger systems need to be 
performed to test this procedure and its results. 

Initially, we attempted to set C and a to zero, thus 
fitting functions to the Ornstein-Zernicke form. We set 
a to zero since it was always very small in the previously 
attempted fitting procedure. If this form is correct, then 
one could ideally find S4{k,0;t) by fitting Si{k^q\t) for 
different angles between k and q under the condition that 
one obtains consistent results. We did not obtain consis- 
tent results for S/i{k,Q\t) with this procedure and also 
found that we needed to fix the value of Si{k,Q\t) to 
obtain values of less than 50. Therefore, to obtain 
an estimate for Si[k,Q]t), we choose to fit Io{k,q;t) for 
q < 1.5 to an Ornstein-Zernicke form and then set the 
value of S4{k,0;t) = Io{k,0;t) where Io{k,0,t) is ob- 
tained from the fits. Note that this is consistent with 
our assumption that the limit limg-^g S^^k, q; t) does not 
depend on the angle between k and q. 

If the glass transition is governed by a growing dy- 
namic length scale, then it is expected that for small 
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FIG. 11: The four smallest wave vectors allowed 
due to periodic boundary conditions of the projection 
Io{k,q;Ta)/Io{k,0;Ta) versus q^iso for T = 0.45, 0.47, 0.5, 
0.55, 0.6, and 0.8. The solid line is the scaling function 
1/(1 + x^). The inset shows all calculated wave vectors less 
than 5. 



enough q that S4{k,q;Ta)/S4{k,0;Ta) versus should 
be described by a universal function F(q^) that is inde- 
pendent of temperature Q • To check if this scaling holds 
for /o(A:, Tfj), we plotted Io{k,q;Ta)/Io{k,0;Ta) versus 
q^iso where lo{k,0;t) and ^iso are obtained from the fits 
described above and the scaling function 1/(1+0;^), which 
is shown in Fig. 1111 The subscript iso in S,iso emphasizes 
that this correlation length was obtained from the ori- 
entational average Io{k, q; Tq) of the four-point structure 
factor 5*4 (fc, (f; Tq). It appears that this scaling holds well 
for the small q values, but we will again caution that 
simulations of larger systems need to be performed to 
verify this observation. Shown as the inset to the fig- 
ure is I(){k, q]Ta) / Io{k,0]Ta) versus q(,iso for wave vec- 
tors with a magnitude less than five, and the deviation 
from the scaling behavior is obvious for the larger wave 
vectors. The correlation length obtained from I^ik, q; Tq) 
is on the order of a particle diameter at the larger tem- 
peratures, but grows to about five particle diameters at 
T ~ 0.45. This growth of the correlation length is con- 
sistent with recent results of Berthier and Jack Note, 
however, that at the lowest temperature S,iso is compara- 
ble to the half the length of the simulation cell, which is 
the largest length we expect to be able to extract from 
the simulation without finite size effects. 

With the values of 5*4 (/c, 0; t„) fixed using the fits from 
Io{k, 0; Tq.), we fit S^^k, q; Tq.) where the angle 9 between 
k and q are 0, 45, and 90 degrees to an Ornstein-Zernicke 
form where only the correlation length is allowed to vary. 
We show Si^k, q; Ta)l Si{k, 0; Tq.) versus g^g, where £,e de- 
pends on the angle 9 between k and q, for T = 0.8, 0.6, 
0.55, 0.5, 0.47, and 0.45 in Fig. [TH Only wave vectors 
with a magnitude less than 1.5 arc shown, which corre- 
sponds to the four smallest wave vectors allowed due to 




FIG. 12: (a) Si{k,q;t) / SA{k,Q;t) versus g^e for 61 = 0°, 45°, 
and 90° where 9 is the angle between k and q calculated for 
the temperatures T = 0.8, 0.6, 0.55, 0.5, 0.47, and 0.45. The 
inset shows the length scales obtained from the different fits 
{6 — 0°, triangles; 9 = 45°, squares; 9 = 90°, circles; (,i3o, 
X's) . (b) Dynamic correlation lengths versus the a relaxation 
time for 9 — 0° (triangles), 45° (squares), and 90° (circles). 
The solid lines are fits to ^ ~ rj. The X's and the dashed 
line corresponds to ^iso obtained from the fits of Io{k,q;Ta). 



periodic boundary conditions at each temperature and 
angle. The overlap is very good for the 18 functions 
shown, and shown in the inset to Fig. [12] are the corre- 
lation lengths. They depend on the angle between k and 
q, and the correlation lengths are largest for 9 = 90° and 
smallest for 9 = 0°. Again, we observe that for 9 = 90°, 
the correlation lengths are larger than half the simula- 
tion cell for r = 0.5 (where ^go w 4.5) and lower. This 
strongly suggests that already at T = 0.5 simulations of 
larger systems are needed in order to verify the present 
results. 

In previous studies it has been found that the correla- 
tion length is related to the a relaxation time according 
to a power law. S, ^ t2 @, SHI]- Recently, this behavior 
was rationalized by the inhomogencous mode-coupling 
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theory We fitted the the correlation lengths to a 
power law of the form ar^^ and obtained values ranging 
from 7 0.22 ± 0.01 for 6i = 0° and 7 = 0.18 ± 0.01 for 
e = 90°, Fig. O Also shown in Fig. [D is £,iso obtained 
from Io{k,q;t); in this case we found 70 = 0.21 ± 0.01, 
which is very close to the previously reported value of 
0.22, [23. Using this analysis, we find that the dynamic 
correlation length is not only different for different angles 
between k and q, but they also grow at a different rate 
as the temperature is lowered and the a relaxation time 
increases. The range of correlations for particles moving 
in the same direction are longer than for particles mov- 
ing in different directions, but it increases slower with 
decreasing temperature. 

Another scaling prediction is that S4{k,0;Ta) ~ t^. 
To test this prediction we fit /o(fc, 0; Tq) to the form ar^. 
In this way we obtain A w 0.37, which is again very close 
to the value of 0.4 reported in Ref. [2^. These values are 
slightly smaller than the recent inhomogeneous mode- 
couphng theory prediction of A = 0.5 01 ■ 

V. CONCLUSIONS 

There have been many studies looking for a growing 
length scale that accompanies the drastic slowing down 
of the dynamics in supercooled and glass forming liquids. 
Recently, one such possibility was examined by Biroli et. 
al |2§| where they associated a growing correlation length 
with a point-to-set correlation function in a model su- 
percooled liquid. Ever since the observation of heteroge- 
neous dynamics in supercooled and glassy systems, it has 
been suggested that a dynamic correlation length may be 
associated with the size of the dynamically heterogeneous 
regions. Since two point correlation functions are inad- 
equate to describe the correlated motion of atoms and 
correlated relaxation of the fluid, four-point correlation 
functions have been developed to examine this cooper- 
ative motion. Normally these correlation functions are 
assumed to be isotropic, or arc isotropic by design. How- 
ever, it has been observed that correlated displacements 
of particles are not isotropic, and thus it is not surprising 
that the four-point correlation functions might also not 
be isotropic. 

In this work we examined the anisotropy of a four-point 
correlation function. We found that for distances com- 



parable to the nearest neighbor distance the anisotropy 
initially increases upon supercooling the liquid, but then 
seems to saturate or even decrease at the lowest temper- 
atures. Furthermore, the time scale that this anisotropy 
is a maximum for nearest neighbor distances is around 
the a relaxation time at higher temperatures, but then it 
increases slower with decreasing temperature than the a 
relaxation time and roughly follows the time correspond- 
ing to the peak position of the non-Gaussian parameter 

For larger distances, we also found anisotropy of the 
four-point correlation function. We studied the time 
dependence of this longer ranged anisotropy and found 
that the time at which it is the largest also approxi- 
mately follows Tng in the supercooled liquid. The longer 
range anisotropy introduces a challenge in determining 
the growing dynamic length scale ^ in glass forming sys- 
tems. This difficulty is compounded by the relatively 
small system sizes usually employed in simulational stud- 
ies of the glass transition. We developed a procedure to 
extract effective dynamic length scales, but larger sys- 
tem sizes need to be simulated to verify our results. Our 
procedure suggests that the dynamic correlation length 
is different depending on the relative direction of mo- 
tion of two particles within the fluid. Furthermore, this 
anisotropic length scale also increases at a different rate 
with decreasing temperature. 

We hope that our present work will stimulate future 
research in two different directions. First, we advocate 
the need to study larger systems and toperform serious 
finite-size analysis of the results [10, [3l|, [12 . In par- 
ticular, wc expect that in the small q limit four-point 
structure factor 6*4 (A; , q; t) is isotropic and we thus we ex- 
pect its anisotropic component hik^ q; t) to vanish in the 
small (f limit. These expectations should be confirmed by 
simulations of larger systems. Second, wc hope that this 
work will stimulate a development of a theoretical model 
that describes the anisotropy of four-point dynamic cor- 
relations. 
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